# Adaptation preferences and perceptions of responsibility in climate-vulnerable regions: #
# Results from the first cross-national survey of small island states and territories #

# Analysis and Visualization Replication Code

# Preamble #### 
# Load packages #
library(tidyverse)
library(grid)
library(gridExtra)
library(ggtext)
library(DeclareDesign)
library(scales)
library(patchwork)
library(foreign)
library(dplyr)

# Load colors #
darkergreen <- "#197054"
logogreen <- "#1F9A70"
marigold <- "#E89E22"
darkestgreen <- "#104132"
darkermarigold <- "#B57B1D"
gray <- "#959595"



# Load data #
sids.rep <- read.csv("sids_climate_survey.csv")

# Generate list unweighted and low sample size countries for color coding #
unweighted <- sids.rep %>%
  group_by(country, Alpha.2.code) %>% 
  filter(all(weights_AgeGender == 1))

unweighted_country_list <- unique(unweighted$country)
unweighted_countrycode_list <-unique(unweighted$Alpha.2.code)

sampsizes <- sids.rep %>%
  group_by(country, Alpha.2.code) %>%
  summarise(count=n())

samp_below_fifty <- sampsizes$country[which(sampsizes$count <50)]
samp_below_fifty_countrycode <- sampsizes$Alpha.2.code[which(sampsizes$count <50)]

weight_uncertainty_in_blue <- c(unweighted_countrycode_list, samp_below_fifty_countrycode)

#### Open-Ended Tallies ####

## Note this uses the dataset sids_anyopenended_unweightedsample.csv, which includes all 25,931 individuals who completed at least the first block of the survey, not just the ~14,710 who completed + could be weighted

future_openended <- read.csv("sids_anyopenended_unweightedsample.csv") 

future_openended[c('future_open_good', 'future_open_bad')] <- sapply(future_openended[c('future_open_good', 'future_open_bad')], function(x) tolower(x))

all_words <- c( "climate", "global warming", "sea level", "weather", "drought", "flood",
                "heat", "wildfire", "hurricane", "storm", "environment", "air quality", "water",
                "temperature", "greenhouse", "carbon", "emissions", "pollution", "sustainability",
                "energy", "solar", "wind", "renewable", "biodiversity", "natural disaster",
                # Spanish
                "clima", "calentamiento global", "nivel del mar", "tiempo meteorológico", "sequía", "inundación",
                "calor", "incendio forestal", "huracán", "tormenta", "medio ambiente", "calidad del aire", "agua",
                "temperatura", "efecto invernadero", "carbono", "emisiones", "contaminación", "sostenibilidad",
                "energía", "solar", "eólica", "renovables", "biodiversidad", "desastre natural",
                # French
                "climat", "réchauffement climatique", "niveau de la mer", "météo", "sécheresse", "inondation",
                "chaleur", "incendie de forêt", "ouragan", "tempête", "environnement", "qualité de l’air", " eau ",
                "température", "effet de serre", "carbone", "émissions", "pollution", "durabilité",
                "énergie", "solaire", "éolien", "renouvelables", "biodiversité", "catastrophe naturelle",
                # Portuguese
                "clima", "aquecimento global", "nível do mar", "tempo", "seca", "inundação",
                "calor", "incêndio florestal", "furacão", "tempestade", "ambiente", "qualidade do ar", "água",
                "temperatura", "efeito de estufa", "carbono", "emissões", "poluição", "sustentabilidade",
                "energia", "solar", "eólica", "renovável", "biodiversidade", "desastre natural",
                # Dutch
                "klimaat", "opwarming van de aarde", "zeespiegel", "weer", "droogte", "overstroming",
                "hitte", "bosbrand", "orkaan", "storm", "milieu", "luchtkwaliteit", "water",
                "temperatuur", "broeikasgassen", "koolstof", "uitstoot", "vervuiling", "duurzaamheid",
                "energie", "zonne-energie", "windenergie", "hernieuwbaar", "biodiversiteit", "natuurramp")

future_openended$future_good_climate <- ifelse(grepl(paste(all_words, collapse = "|"), future_openended$future_open_good), 1,0)
future_openended$future_bad_climate <- ifelse(grepl(paste(all_words, collapse = "|"), future_openended$future_open_bad), 1,0)

sum(future_openended$future_open_good!=future_openended$future_open_bad) # num of people who answered question

table(future_openended$future_good_climate) #744
table(future_openended$future_bad_climate) #1484

(744+1484)/19312


# redoing on smaller sample that is focus of analysis (because weights can be generated for those 14,710)

future_openended_check <- sids.rep[,colnames(sids.rep)%in%c('future_open_good', 'future_open_bad', "Alpha.2.code")] 
future_openended_check[c('future_open_good', 'future_open_bad')] <- sapply(future_openended_check[c('future_open_good', 'future_open_bad')], function(x) tolower(x))

future_openended_check$future_good_climate <- ifelse(grepl(paste(all_words, collapse = "|"), future_openended_check$future_open_good), 1,0)
future_openended_check$future_bad_climate <- ifelse(grepl(paste(all_words, collapse = "|"), future_openended_check$future_open_bad), 1,0)

sum(future_openended_check$future_open_good!=future_openended_check$future_open_bad) # num of people who answered question

table(future_openended_check$future_good_climate) #465
table(future_openended_check$future_bad_climate) #954

(465+954)/14710 # alternatively, a higher fraction (465+954)/11398, 12.4% of people who answered question


#### Figures ####

###### Fig. 1. Perceptions of vulnerability to extreme weather (map) ####

## See Adaptation_Map_Figure_Replication.R 

##### Fig. 2. Perceived impact of climate change on select groups and Relationship between objective risk indicators and subjective risk perceptions ####
# Left panel: Climate Impacts 
climate_impacts_df <- sids.rep %>% 
  select(contains("climate_impacts"), country, Alpha.2.code, weights_AgeGender) %>%
  filter(Alpha.2.code != "US-HI")

# Recode to ordinal
climate_impacts_df$climate_impacts_ord_1 <- NA
climate_impacts_df$climate_impacts_ord_1[which(climate_impacts_df$climate_impacts_1 == "Only a little")] <- 2
climate_impacts_df$climate_impacts_ord_1[which(climate_impacts_df$climate_impacts_1 == "A moderate amount")] <- 3
climate_impacts_df$climate_impacts_ord_1[which(climate_impacts_df$climate_impacts_1 == "A great deal")] <- 4
climate_impacts_df$climate_impacts_ord_1[which(climate_impacts_df$climate_impacts_1 == "Not at all")] <- 1

climate_impacts_df$climate_impacts_ord_2 <- NA
climate_impacts_df$climate_impacts_ord_2[which(climate_impacts_df$climate_impacts_2 == "Only a little")] <- 2
climate_impacts_df$climate_impacts_ord_2[which(climate_impacts_df$climate_impacts_2 == "A moderate amount")] <- 3
climate_impacts_df$climate_impacts_ord_2[which(climate_impacts_df$climate_impacts_2 == "A great deal")] <- 4
climate_impacts_df$climate_impacts_ord_2[which(climate_impacts_df$climate_impacts_2 == "Not at all")] <- 1

climate_impacts_df$climate_impacts_ord_3 <- NA
climate_impacts_df$climate_impacts_ord_3[which(climate_impacts_df$climate_impacts_3 == "Only a little")] <- 2
climate_impacts_df$climate_impacts_ord_3[which(climate_impacts_df$climate_impacts_3 == "A moderate amount")] <- 3
climate_impacts_df$climate_impacts_ord_3[which(climate_impacts_df$climate_impacts_3 == "A great deal")] <- 4
climate_impacts_df$climate_impacts_ord_3[which(climate_impacts_df$climate_impacts_3 == "Not at all")] <- 1

climate_impacts_df$climate_impacts_ord_4 <- NA
climate_impacts_df$climate_impacts_ord_4[which(climate_impacts_df$climate_impacts_4 == "Only a little")] <- 2
climate_impacts_df$climate_impacts_ord_4[which(climate_impacts_df$climate_impacts_4 == "A moderate amount")] <- 3
climate_impacts_df$climate_impacts_ord_4[which(climate_impacts_df$climate_impacts_4 == "A great deal")] <- 4
climate_impacts_df$climate_impacts_ord_4[which(climate_impacts_df$climate_impacts_4 == "Not at all")] <- 1

climate_impacts_df$climate_impacts_ord_5 <- NA
climate_impacts_df$climate_impacts_ord_5[which(climate_impacts_df$climate_impacts_5 == "Only a little")] <- 2
climate_impacts_df$climate_impacts_ord_5[which(climate_impacts_df$climate_impacts_5 == "A moderate amount")] <- 3
climate_impacts_df$climate_impacts_ord_5[which(climate_impacts_df$climate_impacts_5 == "A great deal")] <- 4
climate_impacts_df$climate_impacts_ord_5[which(climate_impacts_df$climate_impacts_5 == "Not at all")] <- 1


climate_impacts_avg<-climate_impacts_df %>%
  select(contains("climate_impacts_ord"), Alpha.2.code,  weights_AgeGender) %>%
  gather(var,val,-c(Alpha.2.code,  weights_AgeGender))%>%
  dplyr::filter(val!="") %>%
  mutate(var=as.factor(var))%>%
  group_by(Alpha.2.code, var)%>%
  dplyr::summarise(wm = weighted.mean(val, w = weights_AgeGender))

climate_impacts_avg$var <- ordered(climate_impacts_avg$var, levels=c("climate_impacts_ord_1","climate_impacts_ord_2","climate_impacts_ord_3","climate_impacts_ord_4", "climate_impacts_ord_5"),
                                   labels=c("You personally",  "Your community", "Other communities in your country",   "Developing countries",  "Rich countries"))

climate_impacts_avg <- climate_impacts_avg%>%
  mutate(Alpha.2.code.label = paste("<span style = 'color: ",
                                    ifelse (Alpha.2.code %in% weight_uncertainty_in_blue, "blue", "gray"),
                                    ";'>",
                                    Alpha.2.code,
                                    "</span>", sep = ""),
         Alpha.2.code.label = fct_reorder(Alpha.2.code.label, as.character(Alpha.2.code))) 

ci <- ggplot(climate_impacts_avg, aes(x = var, y = wm, group = var)) +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) +
  geom_text(aes(label = Alpha.2.code), size = 3, position = position_jitter(width = .2, seed = 123), alpha = 0.75, color = gray) +
  scale_x_discrete(labels = label_wrap(20)) +
  scale_y_continuous(name = "Perceived Impact Scale", limits = c(1, 4), breaks = c(1, 2, 3, 4)) +
  labs(x = "", y = "Perceived Impact Scale", title = "CLIMATE CHANGE") +
  theme(
    panel.grid.minor.x = element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(),
    axis.line = element_line(colour = "black"),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_blank(),
    plot.title = element_text(vjust = -2)
  )

# Add in points for the US #
# CCAM Beliefs & Attitudes, Fall 2023 #
# https://osf.io/jw79p/

ccam_climate_harm <- read.csv("ccam_oct_2023.csv")

ccam_climate_harm <- ccam_US_oct_2023 %>% 
  select(harm_personally, harm_US, harm_dev_countries, weight_wave) %>% 
  mutate(Alpha.2.code = "US")

# Recode to ordinal
ccam_climate_harm$harm_personally_ord <- NA
ccam_climate_harm$harm_personally_ord[which(ccam_climate_harm$harm_personally == "Only a little")] <- 2
ccam_climate_harm$harm_personally_ord[which(ccam_climate_harm$harm_personally == "A moderate amount")] <- 3
ccam_climate_harm$harm_personally_ord[which(ccam_climate_harm$harm_personally == "A great deal")] <- 4
ccam_climate_harm$harm_personally_ord[which(ccam_climate_harm$harm_personally == "Not at all")] <- 1

ccam_climate_harm$harm_US_ord <- NA
ccam_climate_harm$harm_US_ord[which(ccam_climate_harm$harm_US == "Only a little")] <- 2
ccam_climate_harm$harm_US_ord[which(ccam_climate_harm$harm_US == "A moderate amount")] <- 3
ccam_climate_harm$harm_US_ord[which(ccam_climate_harm$harm_US == "A great deal")] <- 4
ccam_climate_harm$harm_US_ord[which(ccam_climate_harm$harm_US == "Not at all")] <- 1

ccam_climate_harm$harm_dev_countries_ord <- NA
ccam_climate_harm$harm_dev_countries_ord[which(ccam_climate_harm$harm_dev_countries == "Only a little")] <- 2
ccam_climate_harm$harm_dev_countries_ord[which(ccam_climate_harm$harm_dev_countries == "A moderate amount")] <- 3
ccam_climate_harm$harm_dev_countries_ord[which(ccam_climate_harm$harm_dev_countries == "A great deal")] <- 4
ccam_climate_harm$harm_dev_countries_ord[which(ccam_climate_harm$harm_dev_countries == "Not at all")] <- 1

ccam_climate_harm_avg <- ccam_climate_harm %>%
  select(contains("_ord"), Alpha.2.code, weight_wave) %>%
  gather(var, val, -c(Alpha.2.code, weight_wave)) %>%
  filter(!is.na(val)) %>%
  mutate(var = as.factor(var)) %>%
  group_by(Alpha.2.code, var) %>%
  dplyr::summarise(
    wm = weighted.mean(val, w = weight_wave),
    se = sqrt(sum(weight_wave * (val - weighted.mean(val, w = weight_wave))^2) / ((sum(weight_wave) - 1) * sum(weight_wave))),
    CI_lower = wm - qt(0.975, df = n() - 1) * se,  # 95% confidence interval lower bound
    CI_upper = wm + qt(0.975, df = n() - 1) * se,  # 95% confidence interval upper bound
    .groups = 'drop'
  )

ccam_climate_harm_avg$var <- ordered(ccam_climate_harm_avg$var, levels=c("harm_personally_ord","harm_US_ord","harm_dev_countries_ord"),
                                     labels=c("You personally",  "Other communities in your country",  "Developing countries"))

# Data from the CCAM report
# https://climatecommunication.yale.edu/wp-content/uploads/2024/01/climate-change-american-mind-beliefs-attitudes-fall-2023.pdf
gw_harm_US <- data.frame(
  response = c("A great deal", "A moderate amount", "Only a little", "Not at all", "Don't know"),
  response_ord = c(4,3,2,1,NA),
  gw_harm_personal = c(16, 29, 24, 22, 10),
  gw_harm_community = c(19, 31, 20, 19,11),
  gw_harm_US = c(30, 29, 16, 15, 10),
  gw_harm_developing =c(45, 21, 10, 12, 12),
  stringsAsFactors = FALSE
)

gw_harm_US <- gw_harm_US %>%
  mutate(gw_harm_personal_num = (gw_harm_personal/100)*1033,
         gw_harm_community_num = (gw_harm_community/100)*1033,
         gw_harm_US_num = (gw_harm_US/100)*1033,
         gw_harm_developing_num = (gw_harm_developing/100)*1033)

# Reshape the data to long format, calculate the mean, and return to wide format
gw_harm_US_avg <- gw_harm_US %>%
  select(contains("num"), response_ord, response) %>%
  gather(var,val, -response, -response_ord) %>%
  filter(!is.na(response_ord)) %>%
  group_by(var) %>%
  summarise(
    m = sum(response_ord * val, na.rm = TRUE) / sum(val, na.rm = TRUE)
  )%>%
  ungroup()

gw_harm_US_avg$var <- ordered(gw_harm_US_avg$var, levels=c("gw_harm_personal_num","gw_harm_community_num", "gw_harm_US_num", "gw_harm_developing_num" ),
                              labels=c("You personally",  "Your community", "Other communities in your country",   "Developing countries"))

gw_harm_US_avg <- gw_harm_US_avg %>% filter(var == "Your community")


ci_us_forso <- ggplot(climate_impacts_avg, aes(x = var, y = wm, group = var)) +
  geom_text(aes(label = Alpha.2.code), size = 3, position = position_jitter(width = .2, seed = 123), alpha = 0.75, color = gray) +
  geom_boxplot(alpha = 0, outlier.shape = NA) +
  geom_point(data = ccam_climate_harm_avg, aes(x = var, y = wm), size = 3,  alpha = 1, color = "red") +
  geom_point(data = gw_harm_US_avg, aes(x = var, y = m), size = 3,  alpha = 1, color = "red") + # Impure "Your community" from report data
  geom_errorbar(data = ccam_climate_harm_avg, aes(ymin = CI_lower, ymax = CI_upper), width = 0.06, size = 0.6, color = "red") +
  scale_x_discrete(labels = label_wrap(20)) +
  scale_y_continuous(name = "Perceived Impact Scale", limits = c(1, 4), breaks = c(1, 2, 3, 4), labels = c("Not at all", "Only a little", "A moderate\n amount", "A great deal")) +
  labs(x = "", y = "Perceived Impact Scale", title = "") +
  theme(
    panel.grid.minor.x = element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(),
    axis.line = element_line(colour = "black"),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text=element_text(size=12),
    axis.title=element_text(size=14,face="bold")
  )

ci_us_forso


# Right Panel Subjectve v Objective Risk with World Bank indicators
slr_worry_df <- sids.rep %>% 
  select(slr_worry, country, Alpha.2.code, weights_AgeGender) 

slr_worry_samp<-slr_worry_df %>%

  dplyr::filter(slr_worry!="") %>%
  dplyr::filter(country!="") %>%
  group_by(Alpha.2.code)%>%
  dplyr::summarise(n=sum(weights_AgeGender, na.rm = TRUE))

slr_worry_worse <- slr_worry_df %>% filter(slr_worry == "Not at all worried" |slr_worry == "Not very worried")  %>%
  dplyr::filter(country!="") %>%
  group_by(Alpha.2.code) %>%
  dplyr::summarise(worse_sum_for_order=sum(weights_AgeGender, na.rm = TRUE))

slr_worry_freq<-slr_worry_df %>%
  dplyr::filter(slr_worry!="") %>%
  dplyr::filter(country!="") %>%
  group_by(Alpha.2.code, slr_worry)%>%
  dplyr::summarise(freq=sum(weights_AgeGender, na.rm = TRUE))  %>%
  dplyr::left_join(slr_worry_samp,by="Alpha.2.code") %>%
  mutate(prop=freq/n) %>%
  dplyr::left_join(slr_worry_worse,by="Alpha.2.code") 



slr_worry_freq$slr_worry <- ordered(slr_worry_freq$slr_worry, levels=c("Very worried","Somewhat worried","Not very worried","Not at all worried"),
                                    labels=c("Very worried","Somewhat worried","Not very worried","Not at all worried"))


slr_worry_freq<-slr_worry_freq%>%
  mutate(per=prop*100,
         worse_prop_for_order = worse_sum_for_order/n)

slr_worry_freq$Alpha.2.code<- reorder(slr_worry_freq$Alpha.2.code, slr_worry_freq$worse_prop_for_order)

slr_worry_freq <- slr_worry_freq%>%
  mutate(Alpha.2.code.label = paste("<span style = 'color: ",
                                    ifelse (Alpha.2.code %in% weight_uncertainty_in_blue, "blue", "black"),
                                    ";'>",
                                    Alpha.2.code,
                                    "</span>", sep = ""))

slr_worry_freq$Alpha.2.code.label<- reorder(slr_worry_freq$Alpha.2.code.label, slr_worry_freq$worse_prop_for_order)

wb_pop_below_5m <- read.csv("API_EN.POP.EL5M.ZS_DS2_en_csv_v2_1739411.csv", sep = ",", skip = 4)
countrynames_crosswalk <- read.csv("countrynames_crosswalk.csv") %>% select(name, alpha.2, alpha.3)
alpha2_list <- unique(sids.rep$Alpha.2.code)
wb_pop_below_5m <- wb_pop_below_5m %>% 
  left_join (countrynames_crosswalk, by = c("Country.Code" = "alpha.3")) %>%
  select(Country.Name, name, Country.Code, alpha.2, X2015) %>% 
  filter(alpha.2 %in% alpha2_list ) %>%
  select(alpha.2, X2015)

slr_worry_freq_worried_only <- slr_worry_freq %>% 
  filter(slr_worry == "Somewhat worried" | slr_worry == "Very worried") %>%
  select(Alpha.2.code, per) %>%
  group_by(Alpha.2.code) %>%
  summarise(percent_worried = sum(per))

slr_worry_freq_worried_only <- slr_worry_freq_worried_only %>%
  left_join(wb_pop_below_5m, by = c("Alpha.2.code" = "alpha.2")) %>%
  filter(!is.na(X2015))

# List of Alpha2codes missing from wb data, mostly territories
# setdiff(alpha2_list, unique(wb_pop_below_5m$alpha.2))
# [1] "AI"     "RE"     "GP"     "YT"     "US-HI"  "MQ"    "HN-IB"  "CO-SAP" "BQ"    

so <- ggplot(slr_worry_freq_worried_only, aes(x = percent_worried, y = X2015)) +
  geom_text(aes(label = Alpha.2.code), size =3) +
  labs(x = "% population worried about sea-level rise", y = "% population living in low-lying elevation", title = "") +
  scale_x_continuous(limits = c(0, 100), breaks = seq(0, 100, 25)) +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 25)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dotted") +
  theme(
    panel.grid.minor.x = element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(),
    axis.line = element_line(colour = "black"),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_blank(),
    plot.title = element_text(vjust = -2,),
    axis.text=element_text(size=12),
    axis.title=element_text(size=14,face="bold")
  )

so

# facet climate impacts and subjective vs objective slr risk/worry
ci_so <- ci_us_forso | so
ci_so



##### Fig. 3. Attribution of responsibility to select countries for causing and for solving climate change ####
# Top panel: Responsible for CAUSING climate change 
responsible_cause_df <- sids.rep %>% 
  select(contains("responsible_cause"), contains("cause_ter"), country, Alpha.2.code, ter, col, weights_AgeGender) %>%
  filter(country != "Belize") %>% # Colonial Power NA 
  filter(Alpha.2.code != "US-HI") 

responsible_cause_df <- responsible_cause_df[-which(responsible_cause_df$col==""),] # individuals who inadvertently removed col metadata from url and thus who may have seen an error-filled survey

# Combine/merge responsible_cause and cause_ter
responsible_cause_df$responsible_all_1 <- as.character(responsible_cause_df$responsible_cause_1)
responsible_cause_df$responsible_all_1[which(responsible_cause_df$ter == "yes")] <- as.character(responsible_cause_df$cause_ter_1[which(responsible_cause_df$ter == "yes")])

responsible_cause_df$responsible_all_2 <- as.character(responsible_cause_df$responsible_cause_2)
responsible_cause_df$responsible_all_2[which(responsible_cause_df$ter == "yes")] <- as.character(responsible_cause_df$cause_ter_2[which(responsible_cause_df$ter == "yes")])

responsible_cause_df$responsible_all_3 <- as.character(responsible_cause_df$responsible_cause_3)
responsible_cause_df$responsible_all_3[which(responsible_cause_df$ter == "yes")] <- as.character(responsible_cause_df$cause_ter_3[which(responsible_cause_df$ter == "yes")])

responsible_cause_df$responsible_all_4 <- as.character(responsible_cause_df$responsible_cause_4)
responsible_cause_df$responsible_all_4[which(responsible_cause_df$ter == "yes")] <- as.character(responsible_cause_df$cause_ter_4[which(responsible_cause_df$ter == "yes")])

responsible_cause_df$responsible_all_5 <- as.character(responsible_cause_df$responsible_cause_5)
responsible_cause_df$responsible_all_5[which(responsible_cause_df$ter == "yes")] <- as.character(responsible_cause_df$cause_ter_5[which(responsible_cause_df$ter == "yes")])

responsible_cause_df$responsible_all_6 <- as.character(responsible_cause_df$responsible_cause_6)
responsible_cause_df$responsible_all_6[which(responsible_cause_df$ter == "yes")] <- as.character(responsible_cause_df$cause_ter_6[which(responsible_cause_df$ter == "yes")])

# Recode to ordinal 
responsible_cause_df$responsible_all_ord_1 <- NA
responsible_cause_df$responsible_all_ord_1[which(responsible_cause_df$responsible_all_1 == "Not very responsible")] <- 2
responsible_cause_df$responsible_all_ord_1[which(responsible_cause_df$responsible_all_1 == "Somewhat resonsible")] <- 3
responsible_cause_df$responsible_all_ord_1[which(responsible_cause_df$responsible_all_1 == "Very responsible")] <- 4
responsible_cause_df$responsible_all_ord_1[which(responsible_cause_df$responsible_all_1 == "Not at all responsible`")] <- 1

responsible_cause_df$responsible_all_ord_2 <- NA
responsible_cause_df$responsible_all_ord_2[which(responsible_cause_df$responsible_all_2 == "Not very responsible")] <- 2
responsible_cause_df$responsible_all_ord_2[which(responsible_cause_df$responsible_all_2 == "Somewhat resonsible")] <- 3
responsible_cause_df$responsible_all_ord_2[which(responsible_cause_df$responsible_all_2 == "Very responsible")] <- 4
responsible_cause_df$responsible_all_ord_2[which(responsible_cause_df$responsible_all_2 == "Not at all responsible`")] <- 1

responsible_cause_df$responsible_all_ord_3 <- NA
responsible_cause_df$responsible_all_ord_3[which(responsible_cause_df$responsible_all_3 == "Not very responsible")] <- 2
responsible_cause_df$responsible_all_ord_3[which(responsible_cause_df$responsible_all_3 == "Somewhat resonsible")] <- 3
responsible_cause_df$responsible_all_ord_3[which(responsible_cause_df$responsible_all_3 == "Very responsible")] <- 4
responsible_cause_df$responsible_all_ord_3[which(responsible_cause_df$responsible_all_3 == "Not at all responsible`")] <- 1

responsible_cause_df$responsible_all_ord_4 <- NA
responsible_cause_df$responsible_all_ord_4[which(responsible_cause_df$responsible_all_4 == "Not very responsible")] <- 2
responsible_cause_df$responsible_all_ord_4[which(responsible_cause_df$responsible_all_4 == "Somewhat resonsible")] <- 3
responsible_cause_df$responsible_all_ord_4[which(responsible_cause_df$responsible_all_4 == "Very responsible")] <- 4
responsible_cause_df$responsible_all_ord_4[which(responsible_cause_df$responsible_all_4 == "Not at all responsible`")] <- 1

responsible_cause_df$responsible_all_ord_5 <- NA
responsible_cause_df$responsible_all_ord_5[which(responsible_cause_df$responsible_all_5 == "Not very responsible")] <- 2
responsible_cause_df$responsible_all_ord_5[which(responsible_cause_df$responsible_all_5 == "Somewhat resonsible")] <- 3
responsible_cause_df$responsible_all_ord_5[which(responsible_cause_df$responsible_all_5 == "Very responsible")] <- 4
responsible_cause_df$responsible_all_ord_5[which(responsible_cause_df$responsible_all_5 == "Not at all responsible`")] <- 1

responsible_cause_df$responsible_all_ord_6 <- NA
responsible_cause_df$responsible_all_ord_6[which(responsible_cause_df$responsible_all_6 == "Not very responsible")] <- 2
responsible_cause_df$responsible_all_ord_6[which(responsible_cause_df$responsible_all_6 == "Somewhat resonsible")] <- 3
responsible_cause_df$responsible_all_ord_6[which(responsible_cause_df$responsible_all_6 == "Very responsible")] <- 4
responsible_cause_df$responsible_all_ord_6[which(responsible_cause_df$responsible_all_6 == "Not at all responsible`")] <- 1

# Countries with Japan as col are actually US
responsible_cause_df$colonialpower <- as.character(responsible_cause_df$col)
responsible_cause_df$colonialpower[which(responsible_cause_df$col == "Japan")] <- "US"


responsible_cause_avg<-responsible_cause_df %>% 
  select(contains("responsible_all_ord"), Alpha.2.code, colonialpower, weights_AgeGender) %>%
  gather(var,val,-c(Alpha.2.code, colonialpower, weights_AgeGender))%>%
  dplyr::filter(val!="") %>%
  mutate(var=as.factor(var))%>%
  group_by(Alpha.2.code, colonialpower, var)%>%
  dplyr::summarise(wm = weighted.mean(val, w = weights_AgeGender))

responsible_cause_avg$var <- ordered(responsible_cause_avg$var, levels=c("responsible_all_ord_1","responsible_all_ord_2","responsible_all_ord_3","responsible_all_ord_4", "responsible_all_ord_5", "responsible_all_ord_6"),
                                     labels=c("US",  "China", "Saudi\nArabia",   "Home\nCountry",  "Regional\nPower",  "Colonial\nPower"))
mean(responsible_cause_avg$wm[responsible_cause_avg$var=="Colonial\nPower"]) ## 3.18
mean(responsible_cause_avg$wm[responsible_cause_avg$var=="US"]) ## 3.45
mean(responsible_cause_avg$wm[responsible_cause_avg$var=="China"]) # 3.45
mean(responsible_cause_avg$wm[responsible_cause_avg$var=="Home\nCountry"]) # 2.6

axis_order <- c("US",  "China", "Saudi\nArabia",   "Home\nCountry",  "Regional\nPower",  "Colonial\nPower")

r1 <- ggplot(responsible_cause_avg, aes(x = var, y = wm, group = var)) +
  geom_text(data = subset(responsible_cause_avg, var != "Colonial\nPower"), aes(label = Alpha.2.code), size = 3, position = position_jitter(width = .2, seed = 123), alpha = 0.5, color = "gray") +
  geom_text(data = subset(responsible_cause_avg, var == "Colonial\nPower"), aes(color = colonialpower, label = Alpha.2.code), size = 3, alpha = 0.5, position = position_jitter(width = .2, seed = 123), key_glyph = "rect") +
  geom_boxplot(alpha = 0, outlier.shape = NA) +
  scale_x_discrete(limits = axis_order) +
  scale_y_continuous(name = "", limits = c(1, 4), breaks = c(1, 2, 3, 4), labels = c("Not at all\nresponsible", "Not very\nresponsible", "Somewhat\nresponsible", "Very\nresponsible")) + #
  labs(title = "Responsibility for causing climate change",x = "", color = "Colonial Power") +
  guides(color = guide_legend(title = NULL)) +
  theme(
    panel.grid.minor.x = element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(),
    axis.line = element_line(colour = "black"),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_blank()#,
    #plot.title = element_text(vjust = -2)
  )

r1

# Bottom panel: Responsible for SOLVING climate change 
#># responsible_solve #
responsible_solve_df <- sids.rep %>% 
  select(contains("responsible_solve"), contains("cause_ter"), country, Alpha.2.code, ter, col, weights_AgeGender) %>%
  filter(country != "Belize") %>% # Colonial Power NA 
  filter(Alpha.2.code != "US-HI")

responsible_solve_df <- responsible_solve_df[-which(responsible_solve_df$col==""),]

# Combine/merge responsible_solve and cause_ter
responsible_solve_df$responsible_all_1 <- as.character(responsible_solve_df$responsible_solve_1)
responsible_solve_df$responsible_all_1[which(responsible_solve_df$ter == "yes")] <- as.character(responsible_solve_df$cause_ter_1[which(responsible_solve_df$ter == "yes")])

responsible_solve_df$responsible_all_2 <- as.character(responsible_solve_df$responsible_solve_2)
responsible_solve_df$responsible_all_2[which(responsible_solve_df$ter == "yes")] <- as.character(responsible_solve_df$cause_ter_2[which(responsible_solve_df$ter == "yes")])

responsible_solve_df$responsible_all_3 <- as.character(responsible_solve_df$responsible_solve_3)
responsible_solve_df$responsible_all_3[which(responsible_solve_df$ter == "yes")] <- as.character(responsible_solve_df$cause_ter_3[which(responsible_solve_df$ter == "yes")])

responsible_solve_df$responsible_all_4 <- as.character(responsible_solve_df$responsible_solve_4)
responsible_solve_df$responsible_all_4[which(responsible_solve_df$ter == "yes")] <- as.character(responsible_solve_df$cause_ter_4[which(responsible_solve_df$ter == "yes")])

responsible_solve_df$responsible_all_5 <- as.character(responsible_solve_df$responsible_solve_5)
responsible_solve_df$responsible_all_5[which(responsible_solve_df$ter == "yes")] <- as.character(responsible_solve_df$cause_ter_5[which(responsible_solve_df$ter == "yes")])

responsible_solve_df$responsible_all_6 <- as.character(responsible_solve_df$responsible_solve_6)
responsible_solve_df$responsible_all_6[which(responsible_solve_df$ter == "yes")] <- as.character(responsible_solve_df$cause_ter_6[which(responsible_solve_df$ter == "yes")])

# Recode to ordinal 
responsible_solve_df$responsible_all_ord_1 <- NA
responsible_solve_df$responsible_all_ord_1[which(responsible_solve_df$responsible_all_1 == "Not very responsible")] <- 2
responsible_solve_df$responsible_all_ord_1[which(responsible_solve_df$responsible_all_1 == "Somewhat resonsible")] <- 3
responsible_solve_df$responsible_all_ord_1[which(responsible_solve_df$responsible_all_1 == "Very responsible")] <- 4
responsible_solve_df$responsible_all_ord_1[which(responsible_solve_df$responsible_all_1 == "Not at all responsible`")] <- 1

responsible_solve_df$responsible_all_ord_2 <- NA
responsible_solve_df$responsible_all_ord_2[which(responsible_solve_df$responsible_all_2 == "Not very responsible")] <- 2
responsible_solve_df$responsible_all_ord_2[which(responsible_solve_df$responsible_all_2 == "Somewhat resonsible")] <- 3
responsible_solve_df$responsible_all_ord_2[which(responsible_solve_df$responsible_all_2 == "Very responsible")] <- 4
responsible_solve_df$responsible_all_ord_2[which(responsible_solve_df$responsible_all_2 == "Not at all responsible`")] <- 1

responsible_solve_df$responsible_all_ord_3 <- NA
responsible_solve_df$responsible_all_ord_3[which(responsible_solve_df$responsible_all_3 == "Not very responsible")] <- 2
responsible_solve_df$responsible_all_ord_3[which(responsible_solve_df$responsible_all_3 == "Somewhat resonsible")] <- 3
responsible_solve_df$responsible_all_ord_3[which(responsible_solve_df$responsible_all_3 == "Very responsible")] <- 4
responsible_solve_df$responsible_all_ord_3[which(responsible_solve_df$responsible_all_3 == "Not at all responsible`")] <- 1

responsible_solve_df$responsible_all_ord_4 <- NA
responsible_solve_df$responsible_all_ord_4[which(responsible_solve_df$responsible_all_4 == "Not very responsible")] <- 2
responsible_solve_df$responsible_all_ord_4[which(responsible_solve_df$responsible_all_4 == "Somewhat resonsible")] <- 3
responsible_solve_df$responsible_all_ord_4[which(responsible_solve_df$responsible_all_4 == "Very responsible")] <- 4
responsible_solve_df$responsible_all_ord_4[which(responsible_solve_df$responsible_all_4 == "Not at all responsible`")] <- 1

responsible_solve_df$responsible_all_ord_5 <- NA
responsible_solve_df$responsible_all_ord_5[which(responsible_solve_df$responsible_all_5 == "Not very responsible")] <- 2
responsible_solve_df$responsible_all_ord_5[which(responsible_solve_df$responsible_all_5 == "Somewhat resonsible")] <- 3
responsible_solve_df$responsible_all_ord_5[which(responsible_solve_df$responsible_all_5 == "Very responsible")] <- 4
responsible_solve_df$responsible_all_ord_5[which(responsible_solve_df$responsible_all_5 == "Not at all responsible`")] <- 1

responsible_solve_df$responsible_all_ord_6 <- NA
responsible_solve_df$responsible_all_ord_6[which(responsible_solve_df$responsible_all_6 == "Not very responsible")] <- 2
responsible_solve_df$responsible_all_ord_6[which(responsible_solve_df$responsible_all_6 == "Somewhat resonsible")] <- 3
responsible_solve_df$responsible_all_ord_6[which(responsible_solve_df$responsible_all_6 == "Very responsible")] <- 4
responsible_solve_df$responsible_all_ord_6[which(responsible_solve_df$responsible_all_6 == "Not at all responsible`")] <- 1

# Countries with Japan as col are actually US
responsible_solve_df$colonialpower <- as.character(responsible_solve_df$col)
responsible_solve_df$colonialpower[which(responsible_solve_df$col == "Japan")] <- "US"


responsible_solve_avg<-responsible_solve_df %>% 
  select(contains("responsible_all_ord"), Alpha.2.code, colonialpower, weights_AgeGender) %>%
  gather(var,val,-c(Alpha.2.code, colonialpower, weights_AgeGender))%>%
  dplyr::filter(val!="") %>%
  mutate(var=as.factor(var))%>%
  group_by(Alpha.2.code, colonialpower, var)%>%
  dplyr::summarise(wm = weighted.mean(val, w = weights_AgeGender))

responsible_solve_avg$var <- ordered(responsible_solve_avg$var, levels=c("responsible_all_ord_1","responsible_all_ord_2","responsible_all_ord_3","responsible_all_ord_4", "responsible_all_ord_5", "responsible_all_ord_6"),
                                     labels=c("US",  "China", "Saudi\nArabia",   "Home\nCountry",  "Regional\nPower",  "Colonial\nPower"))
## average of weighted means 
mean(responsible_solve_avg$wm[responsible_solve_avg$var=="US"]) ## 3.4 
mean(responsible_solve_avg$wm[responsible_solve_avg$var=="China"]) ## 3.3 
mean(responsible_solve_avg$wm[responsible_solve_avg$var=="Colonial\nPower"]) ## 3.2
mean(responsible_solve_avg$wm[responsible_solve_avg$var=="Home\nCountry"],na.rm=TRUE) ## 2.89

r2 <- ggplot(responsible_solve_avg, aes(x = var, y = wm, group = var)) +
  geom_text(data = subset(responsible_solve_avg, var != "Colonial\nPower"), aes(label = Alpha.2.code), size = 3, position = position_jitter(width = .2, seed = 123), alpha = 0.5, color = gray) +
  geom_text(data = subset(responsible_solve_avg, var == "Colonial\nPower"), aes(color = colonialpower, label = Alpha.2.code), size = 3, alpha=0.5, position = position_jitter(width = .2, seed = 123), key_glyph = "rect") +
  geom_boxplot(alpha = 0, outlier.shape = NA) +
  scale_x_discrete(limits = axis_order) +
  scale_y_continuous(name = "", limits = c(1, 4), breaks = c(1, 2, 3, 4),  labels = c("Not at all\nresponsible", "Not very\nresponsible", "Somewhat\nresponsible", "Very\nresponsible")) + # name = "Mean Problem Responsibility"
  labs(title = "Responsibility for solving climate change", x = "",  color = "Colonial Power") +
  guides(color = guide_legend(title = NULL)) +
  theme(
    panel.grid.minor.x = element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(),
    axis.line = element_line(colour = "black"),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_blank()#,
    #plot.title = element_text(vjust = -2)
  )
r2

# Facet/combine to create final figure 
r3 <-  (r1 + theme(legend.position = "none")) / r2 + 
  plot_layout(guides = "collect") & theme(legend.position = "bottom")
r3

y_axis_label <- wrap_elements(grid::textGrob("Mean Responsibility", rot = 90, gp = grid::gpar(fontsize = 14, fontface = "bold")))


# Combine the label and the plots
r3 <- (y_axis_label | r3) + plot_layout(widths = c(1, 20))
r3



##### Fig. 4. Perceived need for climate finance contributions by specific countries ####
sids.rep$emissions_who_mod <- as.character(sids.rep$emissions_who)
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Bengladesh")] <- "Bangladesh"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Chine")] <- "China"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Estados Unidos")] <- "UK"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Grã Bretanha")] <- "UK"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Gran Bretaña")] <- "UK"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Great Britain")] <- "UK"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Groot Brittanië")] <- "UK"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Inde")] <- "India"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Índia")] <- "India"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Le Royaume-Uni")] <- "UK"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Les États Unis")] <- "USA"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Suíça")] <- "Switzerland"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Suisse")] <- "Switzerland"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Suiza")] <- "Switzerland"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Switzlerland")] <- "Switzerland"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "The United Kingdom")] <- "UK"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "The United States")] <- "UK"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "United States")] <- "USA"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Verenigde Staten")] <- "USA"
sids.rep$emissions_who_mod[which(sids.rep$emissions_who_mod == "Zwitserland")] <- "Switzerland"
sids.rep$emissions_who_mod <- as.factor(sids.rep$emissions_who_mod)


sids.rep$contribution_agg <- NA
sids.rep$contribution_agg[which(sids.rep$emissions_version=="both")] <- sids.rep$fund_contribite_both_1[which(sids.rep$emissions_version=="both")] 
sids.rep$contribution_agg[which(sids.rep$emissions_version=="history")] <- sids.rep$fund_contribute_hist_1[which(sids.rep$emissions_version=="history")]
sids.rep$contribution_agg[which(sids.rep$emissions_version=="now")] <- sids.rep$fund_contribute_now_1[which(sids.rep$emissions_version=="now")]
sids.rep$contribution_agg[which(sids.rep$emissions_version=="control")] <- sids.rep$fund_contribute_cont_1[which(sids.rep$emissions_version=="control")]

sids.rep$emissions_version <- factor(sids.rep$emissions_version, levels=c("control", "now", "history", "both"))

finance_exper_pooled1 <- tidy(lm_robust(contribution_agg ~ emissions_version, clusters = Alpha.2.code, data = sids.rep, weights=weights_AgeGender))

finance_exper_USA <- tidy(lm_robust(contribution_agg ~ emissions_version, clusters = Alpha.2.code, data = sids.rep[which(sids.rep$emissions_who_mod=="USA"),], weights=weights_AgeGender))
finance_exper_UK <- tidy(lm_robust(contribution_agg ~ emissions_version, clusters = Alpha.2.code, data = sids.rep[which(sids.rep$emissions_who_mod=="UK"),], weights=weights_AgeGender))
finance_exper_Switzerland <- tidy(lm_robust(contribution_agg ~ emissions_version, clusters = Alpha.2.code, data = sids.rep[which(sids.rep$emissions_who_mod=="Switzerland"),], weights=weights_AgeGender))
finance_exper_India <- tidy(lm_robust(contribution_agg ~ emissions_version, clusters = Alpha.2.code, data = sids.rep[which(sids.rep$emissions_who_mod=="India"),], weights=weights_AgeGender))
finance_exper_China <- tidy(lm_robust(contribution_agg ~ emissions_version, clusters = Alpha.2.code, data = sids.rep[which(sids.rep$emissions_who_mod=="China"),], weights=weights_AgeGender))
finance_exper_Bangladesh <- tidy(lm_robust(contribution_agg ~ emissions_version, clusters = Alpha.2.code, data = sids.rep[which(sids.rep$emissions_who_mod=="Bangladesh"),], weights=weights_AgeGender))

# Figure for results
finance_exper_USA <- finance_exper_USA[2:4,] %>%
  #filter(grepl("adapt_",term))%>%
  mutate(ci.high=estimate+1.64*std.error,
         ci.low=estimate-1.64*std.error)%>%
  ggplot(aes(x=estimate,y=term))+
  ggtitle("USA")+
  geom_pointrange(aes(xmin=ci.low,xmax=ci.high))+
  geom_vline(xintercept=0,lty="dotted")+
  scale_y_discrete(labels = c("emissions_versionnow" = "Current Emissions %",
                              "emissions_versionhistory" = "Cumulative Emissions %",
                              "emissions_versionboth" = "Both %")) +
  xlim(-15,10)+
  labs(x="",y="")+
  theme(legend.position = "none",
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
        axis.text.x = element_text(size=11), 
        axis.title.x = element_text(size=13), 
        axis.line.y = element_line(colour = 'black', size=0.5, linetype='solid'),
        axis.text.y = element_text(size=11), 
        axis.title.y = element_text(size=13))

finance_exper_UK <- finance_exper_UK[2:4,] %>%
  mutate(ci.high=estimate+1.64*std.error,
         ci.low=estimate-1.64*std.error)%>%
  ggplot(aes(x=estimate,y=term))+
  ggtitle("UK")+
  geom_pointrange(aes(xmin=ci.low,xmax=ci.high))+
  geom_vline(xintercept=0,lty="dotted")+
  xlim(-15,10)+
  scale_y_discrete(labels = c("emissions_versionnow" = "Current Emissions %",
                              "emissions_versionhistory" = "Cumulative Emissions %",
                              "emissions_versionboth" = "Both %")) +
  labs(x="",y="")+
  theme(legend.position = "none",
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
        axis.text.x = element_text(size=11), 
        axis.title.x = element_text(size=13), 
        axis.line.y = element_line(colour = 'black', size=0.5, linetype='solid'),
        axis.text.y = element_text(size=11), 
        axis.title.y = element_text(size=13))

finance_exper_Switzerland <- finance_exper_Switzerland[2:4,] %>%
  mutate(ci.high=estimate+1.64*std.error,
         ci.low=estimate-1.64*std.error)%>%
  ggplot(aes(x=estimate,y=term))+
  geom_pointrange(aes(xmin=ci.low,xmax=ci.high))+
  geom_vline(xintercept=0,lty="dotted")+
  xlim(-15,10)+
  ggtitle("Switzerland")+
  scale_y_discrete(labels = c("emissions_versionnow" = "Current Emissions %",
                              "emissions_versionhistory" = "Cumulative Emissions %",
                              "emissions_versionboth" = "Both %")) +
  labs(x="",y="")+
  theme(legend.position = "none",
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
        axis.text.x = element_text(size=11), 
        axis.title.x = element_text(size=13), 
        axis.line.y = element_line(colour = 'black', size=0.5, linetype='solid'),
        axis.text.y = element_text(size=11), 
        axis.title.y = element_text(size=13))

finance_exper_Bangladesh <- finance_exper_Bangladesh[2:4,] %>%
  mutate(ci.high=estimate+1.64*std.error,
         ci.low=estimate-1.64*std.error)%>%
  ggplot(aes(x=estimate,y=term))+
  ggtitle("Bangladesh")+
  geom_pointrange(aes(xmin=ci.low,xmax=ci.high))+
  geom_vline(xintercept=0,lty="dotted")+
  scale_y_discrete(labels = c("emissions_versionnow" = "Current Emissions %",
                              "emissions_versionhistory" = "Cumulative Emissions %",
                              "emissions_versionboth" = "Both %")) +
  labs(x="",y="")+
  xlim(-15,10)+
  theme(legend.position = "none",
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
        axis.text.x = element_text(size=11), 
        axis.title.x = element_text(size=13), 
        axis.line.y = element_line(colour = 'black', size=0.5, linetype='solid'),
        axis.text.y = element_text(size=11), 
        axis.title.y = element_text(size=13))

finance_exper_China <- finance_exper_China[2:4,] %>%
  mutate(ci.high=estimate+1.64*std.error,
         ci.low=estimate-1.64*std.error)%>%
  ggplot(aes(x=estimate,y=term))+
  geom_pointrange(aes(xmin=ci.low,xmax=ci.high))+
  geom_vline(xintercept=0,lty="dotted")+
  ggtitle("China")+
  scale_y_discrete(labels = c("emissions_versionnow" = "Current Emissions %",
                              "emissions_versionhistory" = "Cumulative Emissions %",
                              "emissions_versionboth" = "Both %")) +
  labs(x="",y="")+
  xlim(-15,10)+
  theme(legend.position = "none",
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
        axis.text.x = element_text(size=11), 
        axis.title.x = element_text(size=13), 
        axis.line.y = element_line(colour = 'black', size=0.5, linetype='solid'),
        axis.text.y = element_text(size=11), 
        axis.title.y = element_text(size=13))



finance_exper_India <- finance_exper_India[2:4,] %>%
  mutate(ci.high=estimate+1.64*std.error,
         ci.low=estimate-1.64*std.error)%>%
  ggplot(aes(x=estimate,y=term))+
  geom_pointrange(aes(xmin=ci.low,xmax=ci.high))+
  geom_vline(xintercept=0,lty="dotted")+
  scale_y_discrete(labels = c("emissions_versionnow" = "Current Emissions %",
                              "emissions_versionhistory" = "Cumulative Emissions %",
                              "emissions_versionboth" = "Both %")) +
  ggtitle("India")+
  labs(x="",y="")+
  xlim(-15,10)+
  theme(legend.position = "none",
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
        axis.text.x = element_text(size=11), 
        axis.title.x = element_text(size=13), 
        axis.line.y = element_line(colour = 'black', size=0.5, linetype='solid'),
        axis.text.y = element_text(size=11), 
        axis.title.y = element_text(size=13))



grid.arrange(finance_exper_USA, finance_exper_China, finance_exper_India, finance_exper_UK, finance_exper_Switzerland, finance_exper_Bangladesh)


#### Supplementary Figures ####

##### SI Fig. 2. Perceptions of future community conditions, by country or territory #### 
# AND MAIN TEXT DESCRIPTIVE DATA

future_scale_df <- sids.rep %>% 
  select(future_scale, country, Alpha.2.code, weights_AgeGender) %>% 
  filter(country != "Sao") %>%
  filter(Alpha.2.code != "US-HI")

future_scale_samp<-future_scale_df %>%
  dplyr::filter(future_scale!="") %>%
  dplyr::filter(country!="") %>%
  group_by(Alpha.2.code)%>%
  dplyr::summarise(n=sum(weights_AgeGender, na.rm = TRUE))

future_scale_worse <- future_scale_df %>% filter(future_scale == "A bit worse off" |future_scale == "Much worse off")  %>%
  dplyr::filter(country!="") %>%
  group_by(Alpha.2.code) %>%
  dplyr::summarise(worse_sum_for_order=sum(weights_AgeGender, na.rm = TRUE)) 

future_scale_freq<-future_scale_df %>%
  dplyr::filter(future_scale!="") %>%
  dplyr::filter(country!="") %>%
  group_by(Alpha.2.code, future_scale)%>%
  dplyr::summarise(freq=sum(weights_AgeGender, na.rm = TRUE))  %>%
  dplyr::left_join(future_scale_samp,by="Alpha.2.code") %>%
  mutate(prop=freq/n) %>%
  dplyr::left_join(future_scale_worse,by="Alpha.2.code") 

future_scale_freq$future_scale <- ordered(future_scale_freq$future_scale, levels=c( "Much better off",  "A bit better off", "About the same",   "A bit worse off",  "Much worse off"),
                                          labels=c("Much better off",  "A bit better off", "About the same",   "A bit worse off",  "Much worse off"))


future_scale_freq<-future_scale_freq%>%
  mutate(per=prop*100,
         worse_prop_for_order = worse_sum_for_order/n)


# color x axis countries that are unweighted or sample size below 50 in blue
future_scale_freq <- future_scale_freq%>%
  mutate(Alpha.2.code.label = paste("<span style = 'color: ",
                                    ifelse (Alpha.2.code %in% weight_uncertainty_in_blue, "blue", "black"),
                                    ";'>",
                                    Alpha.2.code,
                                    "</span>", sep = "")) 

future_scale_freq$Alpha.2.code.label<- reorder(future_scale_freq$Alpha.2.code.label, future_scale_freq$worse_prop_for_order)


f<-ggplot(future_scale_freq, aes(x=Alpha.2.code.label, fill=future_scale, y=prop)) +
  geom_bar(position=position_stack(vjust=0.5),stat="identity", width = 0.5) +
  geom_text(position=position_stack(vjust=0.5),aes(label=round(per)), size=2, color="white")+ 
  scale_y_continuous(#breaks=c(0,0.5,1),labels=c("0%","50%","100%"))+
    name="",expand = c(0,0), breaks=c(0,0.5,1), labels=c("0%", "50%", "100%")) +
  scale_x_discrete(name="") +
  scale_fill_manual(values=c( "Much worse off" = darkermarigold, "A bit worse off" = marigold, "About the same" = gray,  "A bit better off" = logogreen, "Much better off" = darkergreen))+
  theme(plot.title = element_text(size = 14, face = "bold"), # controls title font + size
        legend.position="top", legend.justification="left", # moves legend to the top and makes it stacked
        legend.title=element_blank(),
        legend.direction = "horizontal",
        legend.spacing.x = unit(0.1, 'cm'),
        axis.text.x=element_markdown(size=8, colour = "black"), #, colour = "black"
        axis.text.y=element_markdown(size=8)) + #, colour = "black"
  theme( # remove the vertical grid lines
    panel.grid.major.x = element_line(colour="black"),
    panel.grid.minor.x = element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
    axis.line=element_line(colour="black"), # add x and y axis lines
    axis.ticks.y=element_blank(),
    axis.ticks.x=element_blank(),
    plot.title=element_text(vjust=-2),
    plot.margin = margin(0, 0.5, 2, 0.25, "cm"))+
  guides(fill = guide_legend(nrow = 1, reverse=TRUE))+
  coord_flip(clip="off") 

#main text descriptives
future_scale_freq$prop[which(future_scale_freq$Alpha.2.code=="AW"&future_scale_freq$future_scale=="Much worse off")]+
  future_scale_freq$prop[which(future_scale_freq$Alpha.2.code=="AW"&future_scale_freq$future_scale=="A bit worse off")]

future_scale_freq$prop[which(future_scale_freq$Alpha.2.code=="KN"&future_scale_freq$future_scale=="Much worse off")]+
  future_scale_freq$prop[which(future_scale_freq$Alpha.2.code=="KN"&future_scale_freq$future_scale=="A bit worse off")]

future_scale_freq$prop[which(future_scale_freq$Alpha.2.code=="MH"&future_scale_freq$future_scale=="Much worse off")]+
  future_scale_freq$prop[which(future_scale_freq$Alpha.2.code=="MH"&future_scale_freq$future_scale=="A bit worse off")]

future_scale_freq$prop[which(future_scale_freq$Alpha.2.code=="NR"&future_scale_freq$future_scale=="Much worse off")]+
  future_scale_freq$prop[which(future_scale_freq$Alpha.2.code=="NR"&future_scale_freq$future_scale=="A bit worse off")]

future_scale_freq$prop[which(future_scale_freq$Alpha.2.code=="MV"&future_scale_freq$future_scale=="Much better off")]+
  future_scale_freq$prop[which(future_scale_freq$Alpha.2.code=="MV"&future_scale_freq$future_scale=="A bit better off")]

future_scale_freq$prop[which(future_scale_freq$Alpha.2.code=="TV"&future_scale_freq$future_scale=="Much better off")]+
  future_scale_freq$prop[which(future_scale_freq$Alpha.2.code=="TV"&future_scale_freq$future_scale=="A bit better off")]

future_scale_freq$prop[which(future_scale_freq$Alpha.2.code=="KI"&future_scale_freq$future_scale=="Much better off")]+
  future_scale_freq$prop[which(future_scale_freq$Alpha.2.code=="KI"&future_scale_freq$future_scale=="A bit better off")]


##### SI Fig. 3. Belief in climate change (map) ####

# See Adaptation_Map_Replication_Figure.R script


##### SI Fig. 4. Perceived exposure to coastal vulnerability ####
slr_worry_df <- sids.rep %>% 
  select(slr_worry, country, Alpha.2.code, weights_AgeGender) 

slr_worry_samp<-slr_worry_df %>%
  dplyr::filter(slr_worry!="") %>%
  dplyr::filter(country!="") %>%
  group_by(Alpha.2.code)%>%
  dplyr::summarise(n=sum(weights_AgeGender, na.rm = TRUE))

slr_worry_worse <- slr_worry_df %>% filter(slr_worry == "Not at all worried" |slr_worry == "Not very worried")  %>%
  dplyr::filter(country!="") %>%
  group_by(Alpha.2.code) %>%
  dplyr::summarise(worse_sum_for_order=sum(weights_AgeGender, na.rm = TRUE))

slr_worry_freq<-slr_worry_df %>%
  dplyr::filter(slr_worry!="") %>%
  dplyr::filter(country!="") %>%
  group_by(Alpha.2.code, slr_worry)%>%
  dplyr::summarise(freq=sum(weights_AgeGender, na.rm = TRUE))  %>%
  dplyr::left_join(slr_worry_samp,by="Alpha.2.code") %>%
  mutate(prop=freq/n) %>%
  dplyr::left_join(slr_worry_worse,by="Alpha.2.code") 



slr_worry_freq$slr_worry <- ordered(slr_worry_freq$slr_worry, levels=c("Very worried","Somewhat worried","Not very worried","Not at all worried"),
                                    labels=c("Very worried","Somewhat worried","Not very worried","Not at all worried"))


slr_worry_freq<-slr_worry_freq%>%
  mutate(per=prop*100,
         worse_prop_for_order = worse_sum_for_order/n)

slr_worry_freq$Alpha.2.code<- reorder(slr_worry_freq$Alpha.2.code, slr_worry_freq$worse_prop_for_order)

# color x axis countries that are unweighted or sample size below 50 in blue
slr_worry_freq <- slr_worry_freq%>%
  mutate(Alpha.2.code.label = paste("<span style = 'color: ",
                                    ifelse (Alpha.2.code %in% weight_uncertainty_in_blue, "blue", "black"),
                                    ";'>",
                                    Alpha.2.code,
                                    "</span>", sep = ""))

slr_worry_freq$Alpha.2.code.label<- reorder(slr_worry_freq$Alpha.2.code.label, slr_worry_freq$worse_prop_for_order)


w1<-ggplot(slr_worry_freq, aes(x=Alpha.2.code.label, fill=slr_worry, y=prop)) +
  geom_bar(position=position_stack(vjust=0.5),stat="identity", width = 0.5) +
  geom_text(position=position_stack(vjust=0.5),aes(label=round(per)), size=2, color="white")+ 
  scale_y_continuous(#breaks=c(0,0.5,1),labels=c("0%","50%","100%"))+
    name="",expand = c(0,0), breaks=c(0,0.5,1), labels=c("0%", "50%", "100%")) +
  scale_x_discrete(name="") +
  scale_fill_manual(values=c("Very worried"=darkergreen,"Somewhat worried"=logogreen,"Not very worried"=marigold,"Not at all worried"= darkermarigold))+
  theme(plot.title = element_text(size = 14, face = "bold"), # controls title font + size
        legend.position="top", legend.justification="left", # moves legend to the top and makes it stacked
        legend.title=element_blank(),
        legend.direction = "horizontal",
        legend.spacing.x = unit(0.1, 'cm'),
        axis.text.x=element_markdown(size=8, colour = "black"), #, colour = "black"
        axis.text.y=element_markdown(size=8)) + #, colour = "black"
  theme( # remove the vertical grid lines
    panel.grid.major.x = element_line(colour="black"),
    panel.grid.minor.x = element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
    axis.line=element_line(colour="black"), # add x and y axis lines
    axis.ticks.y=element_blank(),
    axis.ticks.x=element_blank(),
    plot.title=element_text(vjust=-2),
    plot.margin = margin(0, 0.5, 2, 0.25, "cm"))+
  guides(fill = guide_legend(nrow = 1, reverse=TRUE))+
  coord_flip(clip="off") 


w1



##### SI Fig. 5. Perceived exposure to water contamination ####
slr_water_df <- sids.rep %>% 
  select(slr_water, country, Alpha.2.code, weights_AgeGender)  %>%
  filter(Alpha.2.code != "US-HI")

slr_water_samp<-slr_water_df %>%
  dplyr::filter(slr_water!="") %>%
  dplyr::filter(country!="") %>%
  group_by(Alpha.2.code)%>%
  dplyr::summarise(n=sum(weights_AgeGender, na.rm = TRUE))

slr_water_worse <- slr_water_df %>% filter(slr_water == "A great deal"  |slr_water == "A moderate amount")  %>%
  dplyr::filter(country!="") %>%
  group_by(Alpha.2.code) %>%
  dplyr::summarise(worse_sum_for_order=sum(weights_AgeGender, na.rm = TRUE))

slr_water_freq<-slr_water_df %>%
  dplyr::filter(slr_water!="") %>%
  dplyr::filter(country!="") %>%
  group_by(Alpha.2.code, slr_water)%>%
  dplyr::summarise(freq=sum(weights_AgeGender, na.rm = TRUE))  %>%
  dplyr::left_join(slr_water_samp,by="Alpha.2.code") %>%
  mutate(prop=freq/n) %>%
  dplyr::left_join(slr_water_worse,by="Alpha.2.code") 


slr_water_freq$slr_water <- ordered(slr_water_freq$slr_water, levels=c("Not at all", "Only a little" , "Don't know", "A moderate amount", "A great deal"))


slr_water_freq<-slr_water_freq%>%
  mutate(per=prop*100,
         worse_prop_for_order = worse_sum_for_order/n)


slr_water_freq <- slr_water_freq%>%
  mutate(Alpha.2.code.label = paste("<span style = 'color: ",
                                    ifelse (Alpha.2.code %in% weight_uncertainty_in_blue, "blue", "black"),
                                    ";'>",
                                    Alpha.2.code,
                                    "</span>", sep = ""))

slr_water_freq$Alpha.2.code.label<- reorder(slr_water_freq$Alpha.2.code.label, slr_water_freq$worse_prop_for_order)


w<-ggplot(slr_water_freq, aes(x=Alpha.2.code.label, fill=slr_water, y=prop)) +
  geom_bar(position=position_stack(vjust=0.5),stat="identity", width = 0.5) +
  geom_text(position=position_stack(vjust=0.5),aes(label=round(per)), size=2, color="white")+ 
  scale_y_continuous(#breaks=c(0,0.5,1),labels=c("0%","50%","100%"))+
    name="",expand = c(0,0), breaks=c(0,0.5,1), labels=c("0%", "50%", "100%")) +
  scale_x_discrete(name="") +
  scale_fill_manual(values=c("A great deal"=darkergreen, "A moderate amount"=logogreen,   "Only a little"= marigold,  "Not at all" = darkermarigold, "Don't know" = gray))+
  theme(plot.title = element_text(size = 14, face = "bold"), # controls title font + size
        legend.position="top", legend.justification="left", # moves legend to the top and makes it stacked
        legend.title=element_blank(),
        legend.direction = "horizontal",
        legend.spacing.x = unit(0.1, 'cm'),
        axis.text.x=element_markdown(size=8, colour = "black"), #, colour = "black"
        axis.text.y=element_markdown(size=8)) + #, colour = "black"
  theme( # remove the vertical grid lines
    panel.grid.major.x = element_line(colour="black"),
    panel.grid.minor.x = element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
    axis.line=element_line(colour="black"), # add x and y axis lines
    axis.ticks.y=element_blank(),
    axis.ticks.x=element_blank(),
    plot.title=element_text(vjust=-2),
    plot.margin = margin(0, 0.5, 2, 0.25, "cm"))+
  guides(fill = guide_legend(nrow = 1, reverse=TRUE))+
  coord_flip(clip="off")

w



##### SI Fig. 6. Perceived impact of sea-level rise on select groups ####
slr_impacts_df <- sids.rep %>% 
  select(contains("slr_impacts"), country, Alpha.2.code, weights_AgeGender) %>%
  filter(Alpha.2.code != "US-HI")

# Recode to ordinal
slr_impacts_df$slr_impacts_ord_1 <- NA
slr_impacts_df$slr_impacts_ord_1[which(slr_impacts_df$slr_impacts_1 == "Only a little")] <- 2
slr_impacts_df$slr_impacts_ord_1[which(slr_impacts_df$slr_impacts_1 == "A moderate amount")] <- 3
slr_impacts_df$slr_impacts_ord_1[which(slr_impacts_df$slr_impacts_1 == "A great deal")] <- 4
slr_impacts_df$slr_impacts_ord_1[which(slr_impacts_df$slr_impacts_1 == "Not at all")] <- 1

slr_impacts_df$slr_impacts_ord_2 <- NA
slr_impacts_df$slr_impacts_ord_2[which(slr_impacts_df$slr_impacts_2 == "Only a little")] <- 2
slr_impacts_df$slr_impacts_ord_2[which(slr_impacts_df$slr_impacts_2 == "A moderate amount")] <- 3
slr_impacts_df$slr_impacts_ord_2[which(slr_impacts_df$slr_impacts_2 == "A great deal")] <- 4
slr_impacts_df$slr_impacts_ord_2[which(slr_impacts_df$slr_impacts_2 == "Not at all")] <- 1

slr_impacts_df$slr_impacts_ord_3 <- NA
slr_impacts_df$slr_impacts_ord_3[which(slr_impacts_df$slr_impacts_3 == "Only a little")] <- 2
slr_impacts_df$slr_impacts_ord_3[which(slr_impacts_df$slr_impacts_3 == "A moderate amount")] <- 3
slr_impacts_df$slr_impacts_ord_3[which(slr_impacts_df$slr_impacts_3 == "A great deal")] <- 4
slr_impacts_df$slr_impacts_ord_3[which(slr_impacts_df$slr_impacts_3 == "Not at all")] <- 1

slr_impacts_df$slr_impacts_ord_4 <- NA
slr_impacts_df$slr_impacts_ord_4[which(slr_impacts_df$slr_impacts_4 == "Only a little")] <- 2
slr_impacts_df$slr_impacts_ord_4[which(slr_impacts_df$slr_impacts_4 == "A moderate amount")] <- 3
slr_impacts_df$slr_impacts_ord_4[which(slr_impacts_df$slr_impacts_4 == "A great deal")] <- 4
slr_impacts_df$slr_impacts_ord_4[which(slr_impacts_df$slr_impacts_4 == "Not at all")] <- 1

slr_impacts_df$slr_impacts_ord_5 <- NA
slr_impacts_df$slr_impacts_ord_5[which(slr_impacts_df$slr_impacts_5 == "Only a little")] <- 2
slr_impacts_df$slr_impacts_ord_5[which(slr_impacts_df$slr_impacts_5 == "A moderate amount")] <- 3
slr_impacts_df$slr_impacts_ord_5[which(slr_impacts_df$slr_impacts_5 == "A great deal")] <- 4
slr_impacts_df$slr_impacts_ord_5[which(slr_impacts_df$slr_impacts_5 == "Not at all")] <- 1


slr_impacts_avg<-slr_impacts_df %>%
  select(contains("slr_impacts_ord"), Alpha.2.code,  weights_AgeGender) %>%
  gather(var,val,-c(Alpha.2.code,  weights_AgeGender))%>%
  dplyr::filter(val!="") %>%
  mutate(var=as.factor(var))%>%
  group_by(Alpha.2.code, var)%>%
  dplyr::summarise(wm = weighted.mean(val, w = weights_AgeGender))


slr_impacts_avg$var <- ordered(slr_impacts_avg$var, levels=c("slr_impacts_ord_1","slr_impacts_ord_2","slr_impacts_ord_3","slr_impacts_ord_4", "slr_impacts_ord_5"),
                               labels=c("You personally",  "Your community", "Other communities in your country", "Developing countries",  "Rich countries"))


# slr impacts US coastal communities #
# using Replication Data for: Mildenberger, M., Sahn, A., Miljanich, C. et al. Unintended consequences of using maps to communicate sea-level rise. Nat Sustain 7, 1018–1026 (2024). https://doi.org/10.1038/s41893-024-01380-0. Mil
CA_slr <- read_csv('SeaLevelRiseSurvey_California_January+16,+2021_13.54.csv') %>%
  mutate(geo = "CA") %>%
  .[-c(1,2),] %>% 
  select(contains("slharm"), geo)%>% 
  select(-c(slharmtime))

FL_slr <- read_csv('SeaLevelRiseSurvey_Florida_January+16,+2021_13.51.csv') %>%
  mutate(geo = "FL") %>%
  .[-c(1,2),] %>% 
  select(contains("slharm"), geo)%>% 
  select(-c(slharmtime))

NJ_slr <- read_csv('SeaLevelRiseSurvey_NewJersey_January+16,+2021_13.50.csv') %>%
  mutate(geo = "NJ") %>%
  .[-c(1,2),] %>% 
  select(contains("slharm"), geo)%>% 
  select(-c(slharmtime))

VA_slr <- read_csv('SeaLevelRiseSurvey_Virginia_January+16,+2021_13.50.csv') %>%
  mutate(geo = "VA") %>%
  .[-c(1,2),] %>% 
  select(contains("slharm"), geo)%>% 
  select(-c(slharmtime))

US_slr <- CA_slr %>%
  bind_rows(FL_slr) %>%
  bind_rows(NJ_slr) %>%
  bind_rows(VA_slr)



# Recode to ordinal
US_slr$slharmpersonal_ord <- NA
US_slr$slharmpersonal_ord[which(US_slr$slharmpersonal == "Only a little")] <- 2
US_slr$slharmpersonal_ord[which(US_slr$slharmpersonal == "A moderate amount")] <- 3
US_slr$slharmpersonal_ord[which(US_slr$slharmpersonal == "A great deal")] <- 4
US_slr$slharmpersonal_ord[which(US_slr$slharmpersonal == "Not at all")] <- 1

US_slr$slharmcommunity_ord <- NA
US_slr$slharmcommunity_ord[which(US_slr$slharmbayarea == "Only a little")] <- 2
US_slr$slharmcommunity_ord[which(US_slr$slharmbayarea == "A moderate amount")] <- 3
US_slr$slharmcommunity_ord[which(US_slr$slharmbayarea == "A great deal")] <- 4
US_slr$slharmcommunity_ord[which(US_slr$slharmbayarea == "Not at all")] <- 1

US_slr$slharmUS_ord <- NA
US_slr$slharmUS_ord[which(US_slr$slharmUS == "Only a little")] <- 2
US_slr$slharmUS_ord[which(US_slr$slharmUS == "A moderate amount")] <- 3
US_slr$slharmUS_ord[which(US_slr$slharmUS == "A great deal")] <- 4
US_slr$slharmUS_ord[which(US_slr$slharmUS == "Not at all")] <- 1

US_slr$slharmdeveloping_ord <- NA
US_slr$slharmdeveloping_ord[which(US_slr$slharmdeveloping == "Only a little")] <- 2
US_slr$slharmdeveloping_ord[which(US_slr$slharmdeveloping == "A moderate amount")] <- 3
US_slr$slharmdeveloping_ord[which(US_slr$slharmdeveloping == "A great deal")] <- 4
US_slr$slharmdeveloping_ord[which(US_slr$slharmdeveloping == "Not at all")] <- 1

US_slr$slharmfuture_ord <- NA
US_slr$slharmfuture_ord[which(US_slr$slharmfuture == "Only a little")] <- 2
US_slr$slharmfuture_ord[which(US_slr$slharmfuture == "A moderate amount")] <- 3
US_slr$slharmfuture_ord[which(US_slr$slharmfuture == "A great deal")] <- 4
US_slr$slharmfuture_ord[which(US_slr$slharmfuture == "Not at all")] <- 1


US_slr_impacts_avg <- US_slr %>%
  select(contains("_ord")) %>%
  select (-slharmfuture_ord) %>% 
  gather(var, val) %>%
  dplyr::filter(val != "") %>%
  mutate(val = as.numeric(val), var = as.factor(var)) %>%
  group_by(var) %>%
  dplyr::summarise(m = mean(val, na.rm = TRUE),
                   se = sd(val, na.rm = TRUE) / sqrt(n()),
                   CI_lower = m - qt(0.975, df = n() - 1) * se,  # 95% confidence interval lower bound
                   CI_upper = m + qt(0.975, df = n() - 1) * se )  # 95% confidence interval upper bound

US_slr_impacts_avg$var <- ordered(US_slr_impacts_avg$var, levels=c("slharmpersonal_ord", "slharmcommunity_ord", "slharmUS_ord", "slharmdeveloping_ord"),
                                  labels=c("You personally",  "Your community", "Other communities in your country", "Developing countries"))



si_us <- ggplot(slr_impacts_avg, aes(x = var, y = wm, group = var)) +
  geom_text(aes(label = Alpha.2.code), size = 3, position = position_jitter(width = .2, seed = 123), alpha = 0.75, color = gray) +
  geom_boxplot(alpha = 0, outlier.shape = NA) +
  geom_point(data = US_slr_impacts_avg, aes(x = var, y = m), size = 3, color = "red") +
  geom_errorbar(data = US_slr_impacts_avg, aes(ymin = CI_lower, ymax = CI_upper, y = m), width = 0.06, size = 0.6, color = "red") +
  scale_x_discrete(labels = label_wrap(20)) +
  scale_y_continuous(name = "", limits = c(1, 4), breaks = c(1, 2, 3, 4), labels = c("Not at all", "Only a little", "A moderate\n amount", "A great deal")) +
  labs(x = "", y = "", title = "SEA-LEVEL RISE") +
  theme(
    panel.grid.minor.x = element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(),
    axis.line = element_line(colour = "black"),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_blank()#,
    #plot.title = element_text(vjust = -2)
  )

si_us

writefig("slr_impacts_SI", 8, 5)
si_us
dev.off()


##### SI Fig. 7. Attribution of responsibility of China vs. select colonial powers ####
responsible_cause_avg_wide<-responsible_cause_avg %>% 
  spread(var, wm)

responsible_cause_avg_wide <- responsible_cause_avg_wide[which(responsible_cause_avg_wide$colonialpower=="France"|
                                                                 responsible_cause_avg_wide$colonialpower=="Spain"|
                                                                 responsible_cause_avg_wide$colonialpower=="United Kingdom"),]
rsubChina <- ggplot(responsible_cause_avg_wide, aes(x=`Colonial\nPower`, y=China, label = Alpha.2.code)) + 
  geom_text(aes(color = colonialpower), size = 3, position=position_jitter(width=.2, seed = 123), key_glyph = "rect")+ 
  labs(x = "Mean Colonial Power Problem Responsibility", y = "Mean China Problem Responsibility", color = "Colonial Power") +
  geom_abline(intercept=0, slope=1, linetype="dashed")+
  scale_x_continuous(breaks = c(3, 4), limits = c(2.75, 4), labels = c("Somewhat\nresponsible", "Very\nresponsible")) +
  scale_y_continuous(breaks = c(3, 4), limits = c(2.75, 4), labels = c("Somewhat\nresponsible", "Very\nresponsible")) +
  guides(color = guide_legend(title = NULL)) +
  theme( # remove the vertical grid lines 
    #panel.grid.major.x = element_line(colour="black"),
    # panel.grid.major.y = element_line(colour = "grey"),
    panel.grid.minor.x = element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
    axis.line=element_line(colour="black"), # add x and y axis lines
    axis.ticks.y=element_blank(),
    axis.ticks.x=element_blank(),
    plot.title=element_text(vjust=-2)) #+

rsubChina

responsible_solve_avg_wide<-responsible_solve_avg %>% 
  spread(var, wm)

responsible_solve_avg_wide <- responsible_solve_avg_wide[which(responsible_solve_avg_wide$colonialpower=="France"|
                                                                 responsible_solve_avg_wide$colonialpower=="Spain"|
                                                                 responsible_solve_avg_wide$colonialpower=="United Kingdom"),]

rsub2China <- ggplot(responsible_solve_avg_wide, aes(x=`Colonial\nPower`, y=China, label = Alpha.2.code)) + 
  geom_text(aes(color = colonialpower), size = 3, position=position_jitter(width=.2, seed = 123), key_glyph = "rect")+ 
  labs(x = "Mean Colonial Power Solution Responsibility", y = "Mean China Solution Responsibility", color = "Colonial Power") +
  geom_abline(intercept=0, slope=1, linetype="dashed")+
  scale_x_continuous(breaks = c(3, 4), limits = c(2.75, 4),  labels = c("Somewhat\nresponsible", "Very\nresponsible")) +
  scale_y_continuous(breaks = c(3, 4), limits = c(2.75, 4),  labels = c("Somewhat\nresponsible", "Very\nresponsible")) +  
  guides(color = guide_legend(title = NULL)) +
  theme(
    panel.grid.minor.x = element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
    axis.line=element_line(colour="black"), # add x and y axis lines
    axis.ticks.y=element_blank(),
    axis.ticks.x=element_blank(),
    plot.title=element_text(vjust=-2)) #+

rsub2China


rsubChinafacet <-  (rsubChina + theme(legend.position = "none")) / 
  plot_spacer() /  rsub2China + 
  plot_layout(guides = "collect", heights = c(1, 0.1, 1)) & theme(legend.position = "bottom") 

rsubChinafacet


##### SI Fig. 8. Attribution of responsibility of US vs. select colonial powers ####
rsubUS <- ggplot(responsible_cause_avg_wide, aes(x=`Colonial\nPower`, y=US, label = Alpha.2.code)) + 
  geom_text(aes(color = colonialpower), size = 3, position=position_jitter(width=.2, seed = 123), key_glyph = "rect")+ 
  labs(x = "Mean Colonial Power Problem Responsibility", y = "Mean US Problem Responsibility", color = "Colonial Power") +
  geom_abline(intercept=0, slope=1, linetype="dashed")+
  scale_x_continuous(breaks = c(3, 4), limits = c(2.75, 4), labels = c("Somewhat\nresponsible", "Very\nresponsible")) +
  scale_y_continuous(breaks = c(3, 4), limits = c(2.75, 4), labels = c("Somewhat\nresponsible", "Very\nresponsible")) +
  guides(color = guide_legend(title = NULL)) +
  theme( # remove the vertical grid lines 
    #panel.grid.major.x = element_line(colour="black"),
    # panel.grid.major.y = element_line(colour = "grey"),
    panel.grid.minor.x = element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
    axis.line=element_line(colour="black"), # add x and y axis lines
    axis.ticks.y=element_blank(),
    axis.ticks.x=element_blank(),
    plot.title=element_text(vjust=-2)) #+

rsubUS


rsub2US <- ggplot(responsible_solve_avg_wide, aes(x=`Colonial\nPower`, y=US, label = Alpha.2.code)) + 
  geom_text(aes(color = colonialpower), size = 3, position=position_jitter(width=.2, seed = 123), key_glyph = "rect")+ 
  labs(x = "Mean Colonial Power Solution Responsibility", y = "Mean US Solution Responsibility", color = "Colonial Power") +
  geom_abline(intercept=0, slope=1, linetype="dashed")+
  scale_x_continuous(breaks = c(3, 4), limits = c(2.75, 4),  labels = c("Somewhat\nresponsible", "Very\nresponsible")) +
  scale_y_continuous(breaks = c(3, 4), limits = c(2.75, 4),  labels = c("Somewhat\nresponsible", "Very\nresponsible")) +  
  guides(color = guide_legend(title = NULL)) +
  theme(
    panel.grid.minor.x = element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
    axis.line=element_line(colour="black"), # add x and y axis lines
    axis.ticks.y=element_blank(),
    axis.ticks.x=element_blank(),
    plot.title=element_text(vjust=-2)) #+

rsub2US


rsubUSfacet <-  (rsubUS + theme(legend.position = "none")) / 
  plot_spacer() /  rsub2US + 
  plot_layout(guides = "collect", heights = c(1, 0.1, 1)) & theme(legend.position = "bottom") 

rsubUSfacet



##### SI Table 3. Difference in means test for attribution of responsibility for causing and addressing climate change ####

## difference in means test for responsible cause question: create long-form dataset where Y=responsibility and X=entity, and regress responsibility on entity. This is a simple difference in means test without controlling for country of residence of across colonial powers. 
responsible_cause_long<-responsible_cause_df%>%
  mutate(id=row_number())%>%
  select(id,weights_AgeGender,contains("responsible_all_ord"))%>%
  pivot_longer(cols=contains("responsible"),names_to="country",values_to="responsible_cause")%>%
  mutate(country=gsub("responsible_all_ord_","",country),
         country=case_when(country=='1'~"US", ## extracted these levels from line 2022 below. 
                           country=='2'~"China",
                           country=='3'~"Saudi Arabia",
                           country=='4'~"Home country",
                           country=='5'~"Regional Power",
                           country=='6'~"Colonial Power")#,
  )
responsible_cause_long$country<-factor(responsible_cause_long$country)
responsible_cause_long$country<-relevel(responsible_cause_long$country,ref="US")


causereg<-lm_robust(responsible_cause~country,data=responsible_cause_long,#weights=weights_AgeGender,
                    clusters=id) ## uses respondent weights, clustered SEs by respondent 
summary(causereg)

# difference of means test ##
responsible_solve_long<-responsible_solve_df%>%
  mutate(id=row_number())%>%
  select(id,weights_AgeGender,contains("responsible_all_ord"))%>%
  pivot_longer(cols=contains("responsible"),names_to="country",values_to="responsible_solve")%>%
  mutate(country=gsub("responsible_all_ord_","",country),
         country=case_when(country=='1'~"US", ## extracted these levels from line 2022 below. 
                           country=='2'~"China",
                           country=='3'~"Saudi Arabia",
                           country=='4'~"Home country",
                           country=='5'~"Regional Power",
                           country=='6'~"Colonial Power")#,
         # responsible_solve=scale(responsible_solve)
  )
responsible_solve_long$country<-factor(responsible_solve_long$country)
responsible_solve_long$country<-relevel(responsible_solve_long$country,ref="US")
solvereg<-lm_robust(responsible_solve~country,data=responsible_solve_long,#weights=weights_AgeGender,
                    clusters=id) ## uses respondent weights, clustered SEs by respondent 
summary(solvereg)



##### SI Table 4. Predicting which respondents view extreme weather as something that will impact their country “only a little” or “not at all ####

future_concern <- sids.rep %>% 
  mutate(id=row_number()) %>% 
  select(id, weights_AgeGender, futureexperience, age, gender,edu, wealthcomp_1,
         indigenous, child, ownhome, ocean_dist, ground_level) %>%
  mutate(
    futureexperience_low = case_when(
      futureexperience == "Only a little" ~ 1,
      futureexperience == "Not at all" ~ 1,
      futureexperience == "A great deal" ~ 0,
      futureexperience == "A moderate amount" ~ 0,
      futureexperience == "Don’t know" ~ 0,
      futureexperience == "" ~ 0,
      TRUE ~ NA_real_  
    )
  )


future_concern$edu_reg <- NA
future_concern$edu_reg[which(future_concern$edu=="Post-graduate")]<- "College or Higher"
future_concern$edu_reg[which(future_concern$edu=="University completed")]<- "College or Higher"
future_concern$edu_reg[which(future_concern$edu=="Post-secondary school other than university")]<- "College or Higher"
future_concern$edu_reg[which(future_concern$edu=="Some university")]<- "College or Higher"
future_concern$edu_reg[which(future_concern$edu=="Secondary school / high school completed")]<- "Highschool"
future_concern$edu_reg[which(future_concern$edu=="Intermediate school or some secondary school / high school ")]<- "Less than High School"
future_concern$edu_reg[which(future_concern$edu=="Primary school completed")]<- "Less than High School"
future_concern$edu_reg[which(future_concern$edu=="Some primary schooling ")]<- "Less than High School"
future_concern$edu_reg[which(future_concern$edu=="Informal schooling ")]<- "Less than High School"
future_concern$edu_reg[which(future_concern$edu=="")]<- NA



library(stargazer)
options(scipen = 999)

low_concern <- lm(futureexperience_low~age + gender +indigenous + child +
                    ownhome+wealthcomp_1+ocean_dist+ground_level, data=future_concern)
stargazer(low_concern)


##### SI Fig. 9. Average desired contributions from each country, for respondents in the control group ####
contribution.df<-sids.rep%>%
  filter(emissions_version=="control",
         !is.na(contribution_agg))%>%
  select(contribution_agg,emissions_who_mod,Alpha.2.code)%>%
  group_by(emissions_who_mod)%>%
  summarise(mean_contrib=mean(contribution_agg),
            se_contrib=sd(contribution_agg)/sqrt(n()))%>%
  mutate(ci.high=mean_contrib+1.96*se_contrib,
         ci.low=mean_contrib-1.96*se_contrib)
contribution.df ### @emma why is there a row here with "" as emissions_who_mod? 
nrow(sids.rep[sids.rep$emissions_who_mod=="",])
## for now eliminate that row

## re-order countries based on mean contribution 
contribution.df$emissions_who_mod<-fct_reorder(factor(contribution.df$emissions_who_mod),contribution.df$mean_contrib)
contrib.fig<-contribution.df%>%
  filter(emissions_who_mod!="")%>%
  ggplot(aes(x=emissions_who_mod,y=mean_contrib))+
  geom_bar(stat="identity",fill=marigold)+
  geom_errorbar(aes(ymin=ci.low,ymax=ci.high),width=.2)+
  labs(x="Contributing country",y="Mean desired contribution ($billions)")+
  theme_bw()

contrib.fig



